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Abstract 

The structure and the electronic properties of a sodium tetrasilicate 
(Na2Si40g) glass were studied by combined Car-Parrinello and classical 
molecular-dynamics simulations. The glass sample was prepared using a 
method recently employed in a study of a silica glass [M. Benoit et al, Euro. 
Phys. J. B 13, 631 (2000)]. First we generated a NS4 glass by classical 
molecular-dynamics and then we took it as the initial configuration of a 
first-principles molecular-dynamics simulation. In the ab initio molecular- 
dynamics simulation, the electronic structure was computed in the frame- 
work of the Kohn-Sham density functional theory within the generalized gra- 
dient approximation using a B-LYP functional. The Car-Parrinello dynamics 
is remarkably stable during the considered trajectory and, as soon as it is 
switched on, some significant structural changes occur. The ab initio descrip- 
tion improves the comparison of the structural characteristics with experi- 
mental data, in particular concerning the Si-0 and Na-0 bond lengths. From 
an electronic point of view, we find that the introduction of the sodium oxide 
in the silica network lowers the band gap and leads to a highly non-localized 
effect on the charges of the network atoms. 
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1. INTRODUCTION 



Over the past three decades, a large number of experimental and theoretical investiga- 
tions has been carried out in order to extract useful information about the structure and the 
physical properties of silicate glasses. These materials have stimulated great interest because 
the advances in understanding them are essential for progress in many scientific fields, such 
as glass and ceramics chemistry, electronics, earth sciences, or the confinement of industrial 
and nuclear wastes. However, in spite of their evident importance, there are still important 
ambiguities in their characterization at a microscopic level. This situation is mainly due to 
the inherent structural disorder present in glasses which complicates the understanding of 
the location of the different elements by comparison with crystalline materials. 

In view of the current experimental findings and theoretical models, it is generally ac- 
cepted that the structure of the most common, widely used and intensely studied glass, pure 
amorphous silica, is composed of SiC>4 tetrahedra which are linked together in a continuous 
random network by strong bonds of bridging oxygens (BO), i.e. oxygen atoms bonded to 
two silicon atoms [|TJ. When a given amount of alkali (Li 2 0, Na 2 or K 2 0), or alkaline earth 
(CaO, MgO) metal oxide is incorporated, the connectivity of this network is decreased by 
the breaking of BO bonds and the formation of non-bridging oxygen atoms (NBO) which are 
bonded to only one silicon atom and weakly bonded to a cation modifier. Consequently the 
local environment around the silicon atoms changes depending on the number of bridging 
oxygen atoms lying within the first shell. This is commonly described by means of the nota- 
tion Q n where n is the number of BO atoms linked to a given Si atom (n=4,3,2,l,0). From a 
technological point of view, the introduction of an alkali metal oxide in silica is accompanied 
by a large reduction of the viscosity, an increase in the thermal expansion coefficient, and 
a change in the refractive index: these physical properties all depend on the bridging to 
non-bridging oxygen ratio in the glass . 

In particular an enormous amount of work has been done on binary sodium silicate 
glasses since almost all commercial silicate glasses and geologic magma are based on this 
family of compositions. They may also serve as useful prototypes to interpret the structural 
and transport properties of more complicated silicate glasses and melts (see for example 
|3|-|| and references therein). 

In addition to physical experiments and to advances in techniques such as EXAFS, 
NMR or XPS, classical molecular-dynamics (MD) simulations have began to provide some 
insight into the three-dimensional structure of these glasses. However one has to notice 
that, even though reasonable results have been reported [|IO|-|i7||, this kind of molecular 
dynamics studies utilizes an interatomic potential generally fixed once and forall throughout 
the simulation. Moreover it is difficult to adjust a classical potential functional form and 
parameters which are able to incorporate both a covalent and an ionic character for the bonds 
existing in the real systems in order to obtain a genuine realistic description. However there 
exists an alternative approach, the ab initio molecular dynamics simulations, where these 
drawbacks are, in principle, eliminated as the interatomic potential is updated at every 
simulation step by first principles calculations. The charge transfers between the NBO and 
the modifier ions can thus be correctly accounted for. Unfortunately one has to pay a price 
for these obvious advantages, namely a small system size and an extremely limited maximum 
length of the trajectory. 
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In this article, we present the microscopic characteristics of a sodium tetrasilicate glass 
Na2Si40g (called hereafter NS4) derived from an ab initio molecular dynamics simulation. 
Numerous experimental studies (NMR [p^-pl 



scopies 



H, XAFS [HIS .. 
and the properties of this glass 

calculations 0-011 or in reverse Monte Carlo (RMC) fits j24] 



neutron diffraction [p2}-|24||, Raman spectro- 
have been carried out in order to elucidate the structure 
The NS4 system has been also analyzed in classical MD 

The attention given to 



this system may be justified by the fact that it can be used as a simple model for more 
complicated aluminosilicate and hydrous silicate glasses. 

The present NS4 glass sample was generated using an approach originally employed for 
obtaining and studying a silica glass sample p9[ in two successive steps. First the liquid 



equilibration, the quench and the glass relaxation are performed classically using a model 
potential having the same functional form as the so-called BKS potential (van Best et al. 
|30| ) for pure silica. In a second step, the resulting glass configuration is relaxed using the 
Car-Parrinello (CP) method JH],^]. As in the silica study, a significant amount of CPU 
time is saved by adopting this strategy since the liquid equilibration, the quench and part 
of the low temperature relaxation are carried out within the framework of classical MD. 

The outline of the paper is as follows: In Sec.[| the methodology is briefly described 
and the details of the particular simulations performed here are given. In Sec.[| the main 
results are presented, including structural properties of the NS4 sample in Sec. |3 A| as well 
as electronic properties in Sec.|3 B|. Section gives the conclusions. 



2. SIMULATION DETAILS 

We started by generating a liquid sample at a temperature around 3500 K containing 
90 atoms (24 silicon atoms, 54 oxygen atoms and 12 sodium atoms) confined in a cubic 
box of edge length 10.81 A corresponding to the experimental NS4 mass density of 2.38 



g/cm 3 f33|. The initial configuration of the NS4 liquid was obtained by melting at 3500 K a 
/3-cristobalite crystal in which SiOj 2 tetrahedra were replaced randomly by Na 2 0^ 2 groups. 

For this system we performed classical molecular dynamics simulations within the mi- 
crocanonical ensemble, with periodic boundary conditions. As mentioned above, we have 
chosen the interaction potential to be a generalization of the BKS potential |30| successfully 



applied for silica and which was initially developed by Kramer et al. |35| in order to 



study zeolites. Recently Horbach et al. |17j have adjusted it for sodium silicate studies by 
introducing an additional short range term for the sodium ion in order to reproduce the 
experimental mass density and the structure factor of the sodium disilicate glass. 

The NS4 liquid was equilibrated during 70 ps and then cooled from 3500 K to ~ 300 K 
at a quench rate of 5 x 10 13 K/s, using a time step of 0.7 fs. The cooling procedure used 
here is identical to the one used in former studies of Si022 [p9j , p4|| and it has been shown 
that it gives access to structural, dynamical and thermal properties in good agreement 
with experiments. The resulting glass was relaxed during another 70 ps and then the final 
classical configuration was used as starting point for the Car-Parrinello ab initio simulation. 
Immediately after switching on the CP dynamics we observed a temperature jump similar 
to the one reported in the silica study which was done following the same approach 
(see Fig. p] where the evolution of the ionic temperature is depicted for the end of the 
classical MD trajectory and the full subsequent CP trajectory). Hence in the CP simulation 
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after a very short relaxation (see the inset of the Fig. [j]), the system heats up and then its 
temperature remains stable around an average value of ~ 410 K (87 K higher than that of 
the classical MD simulation). 

The CP simulation was performed with the ab initio MD code, CPMD developed in 
Stuttgart In the ab initio MD simulations, the electronic structure was treated via 

the Kohn-Sham (KS) formulation of density functional theory |37]] within the generalized 
gradient approximation employing the B-LYP functional [ |3"8"| , |3"9"|] . The KS orbitals were 
expanded in a plane wave basis at the T-point of the supercell up to an energy cutoff of 70 
Ry. Core electrons were not treated explicitly but replaced by atomic pseudopotentials of 
the Goedecker type for silicon jlO] and the Troullier-Martins type for oxygen |OJ], while for 
sodium a Goedecker type semi-core pseudopotential was employed. Preliminary tests have 
been performed to check the adequacy of the choices of the pseudopotentials, exchange and 
correlation functionals and energy cutoff value. The choices of pseudopotentials for Si and 
O and the energy cutoff value led to a good convergence for the structural properties at 
K in cristobalite and a-quartz. The sodium pseudopotential and exchange and correlation 
functionals are justified by bond length estimations carried out on Na2 and Na20 molecules 
and which have been found to be close to experimental values [|42| , f43! or more sophisticated 
theoretical calculations available in the literature ||44| | . 

The CP relaxation was performed for m 6 ps, with a time step of 5.5 a.u. (0.13 fs) and 
a fictitious electronic mass of 800 a.u. used for the integration of the equations of motion. 
The data have been averaged over the last 5 ps of the trajectory. 

The procedure described in the above paragraphs has been carried out a second time 
following exactly the same steps and using the same parameters for both BKS and CP 
simulations, but starting from a different BKS liquid sample. Thus we have generated a 
second NS4 glass sample which was once again remarkably stable during the CP trajectory 
after an identical temperature jump showed up at the beginning of the CP simulation. 
Indeed this second sample heats up from an average temperature of ~ 315 K at the end of 
the BKS relaxation to an average CP temperature of ~ 395 K. Noting that the temperature 
jumps were of the same order for both samples, we may explain this phenomenon as being 
due to the fact that under CP, the samples relaxed towards a potential minimum lower that 
the one calculated during the BKS trajectories and this supplementary potential energy 
has been transformed into ionic kinetic energy. In the next sections, we will discuss the 
structural properties obtained by averaging the data collected for both samples during their 
BKS and CP simulations. 



3. RESULTS AND DISCUSSION 
A. Structure 

In this subsection, we present the structural analysis of both NS4 samples by discussing 
commonly interrelated features like the pair correlation functions, bond angle distributions, 
coordination numbers, Q-species distributions, and bridging to non-bridging oxygen ratio. 
These topics have been investigated for both classical and CP samples in order to understand 
how much the local arrangement of network formers and modifiers is changed when the 
system is treated correctly from a quantum mechanical point of view. 
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We start by remarking that both CP and BKS simulations give the same NBO con- 
centration percentage (22.22%) equal to the sodium-ion concentration which is the value 
predicted theoretically and confirmed by experiments [p0p5|]. We have identified the NBO 



atoms as being either the oxygen atoms having one Si and one Na atom as nearest neigh- 
bours or the oxygen atoms having only one Si at a distance below a Si-0 bond cutoff (taken 
to be 2.0 A, i.e. where the Si-0 pair correlation function becomes zero after the first peak 
- see below) and we have obtained the same results. Analyzing the neighbourhood of the 
NBO atoms, we have systematically found that each of these atoms has two Na atoms at 
almost equal distances, as second and third nearest neighbours. We have also considered the 
nearest-neighbours of each Si atom and we have found that all Si atoms were coordinated 
by 4 oxygens during both runs. 

Further, once the NBO atoms have been identified, we have determined the Q-species 
distributions which are equal to 54.2% Q 4 , 41.6% Q 3 , 4.2% Q 2 for the first sample and to 
58.3% Q 4 , 33.3% Q 3 , 8.3% Q 2 for the second one. We note that the switching on of the 
CP description did not change these percentages. For the NS4 glass, previous 29 Si magic- 
angle spinning nuclear magnetic resonance experiments reported the following results : 55% 
Q 4 , 40% Q 3 , 5% Q 2 by Emerson et al. pf and 50% Q 4 , 48% Q 3 , 2% Q 2 by Maekawa 
et al. [pt]|] . Our simulation results for the first sample compare relatively favorably with 
the experimental data, while, for the second one, there is an excess of Q 4 and Q 2 species. 
Nevertheless, in spite of the limited size of the simulated glass and of the ultrafast quench 
rate, we may point out that there is good overall agreement from the experimental and the 
simulation data for the Q-species distribution. 

From a structural point of view, we note that some significant effects appear as soon as 
the CP dynamics starts and all of them take place on a very short period of time (~ 0.1 - 
0.2 ps) which can be related to the stabilization time of the temperature. We will examine, 
in the next paragraphs, each of these effects occurring when switching from the classical to 
the ab initio description and we will discuss the corresponding average structural properties. 



1. Mean angles and angular distributions 
First, as in the silica glass study f29 |, an important variation occurs in the mean value 



of the inter-tetrahedral angle Si-O-Si, which decreases from 146.1° to 141.8° - see Fig. 0(c) 
which shows the time dependence of the Si-O-Si mean bond-angle during the last 5 ps of the 
classical trajectory and the full CP trajectory. Even if the behaviour was similar for both 
samples, we have chosen to present in Fig. only the data corresponding to one sample. 
Similarly, in the remaining part of this subsection, all the quantities discussed as a function 
of time have been plotted for one sample. All the other characteristics have been averaged 
over both samples. 

In Fig. |^(a) and (b), the time dependencies of the O-Si-0 and Si-O-Na mean bond-angles 
during the same time spans are depicted. We can remark that both descriptions give almost 
identical mean values for the tetrahedral angle O-Si-O, very close to the ideal tetrahedral 
angle 109.47°. Concerning the Si-O-Na angle, we have computed it taking into account that 
each NBO has two Na atoms in its nearest neighbourhood as mentioned above. Consequently 
the depicted values have been taken as the average value of the two corresponding Si-NBO- 
Na linkages. By looking at the time evolution of the mean value of this angle, we note 
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firstly that it presents a slight shift to a lower value during the CP trajectory. But, since 
its fluctuations during the CP simulation are slightly larger compared to the BKS ones, we 
cannot say that the observed shift is a real one without knowing what the origin of these 
increasing fluctuations is. At first glance, one may think that these increasing fluctuations 
would be due to the temperature shift. Hence, in order to check this hypothesis, we have 
classically heated up one of the BKS samples to the CP temperature (i.e. approximately 80 
K higher) and let it relax during 70 ps at this new temperature. Neither the mean value of 
the Si-O-Na angle nor the amplitude of the fluctuations showed any real variation for the 
heated BKS sample. In view of this result, we cannot interpret the increasing of the Si-O-Na 
angle fluctuations between the classical and the ab initio description as being a temperature 
effect. They can be much more probably due to the differences in the nearest-neighbour 
distances observed between the BKS and the CP descriptions (see below). On the other 
hand, these changes may indicate the fact that the ionicity of the Na-NBO bonds is better 
described by the CP approach. The Na-NBO bonds are therefore weaker than in the BKS 
description and this would induce larger fluctuations in the Si-O-Na angle and equally in 
the Na-NBO mean distances, as we will see in the next paragraphs. 

In Fig. we have reported the time-averaged distributions of the tetrahedral O-Si-0 
and inter-tetrahedral Si-O-Si angles as well as that of the Si-O-Na bond angle, evaluated 
for both the BKS and CP trajectories. The O-Si-0 angular distribution (Fig. 0(a)) shows 
in both sharp peak at 109° as for the pure silica and this demonstrates that the 

addition of the sodium oxide does not modify the tetrahedral environment of the Si atoms 
as was already pointed out from previous classical MD calculations [|l"2"l , |T6l and in a RMC fit 
of the NS4 neutron diffraction data p4" |. For the Si-O-Na bond-angle (Fig. 0(b)), we note 
that both distributions have the same average values (i.e. 120° ± 19 in the BKS case and 
120° ± 22 in the CP one), even if their shapes are not really identical. 

From the comparison of the Si-O-Si distributions (Fig. 0(c)), it follows that the CP 



distribution presents a shift towards lower values, already noticed in the silica study |2l] 
and which is of course consistent with the decrease of the mean angle mentioned above. 
Nevertheless both distributions have an asymmetric shape as obtained in the recent RMC 
fits |23| even if their average value (136° ± 16) is slightly lower than those given by the 
present CP and BKS simulations (i.e. 141° ± 15 in the CP case and 145° ± 15 in the BKS 
one). However, in the context of the controversy existing in the literature on the Si-O-Si 
angle distribution for pure silica as well as for other silicates, we note that the average values 
given here are comparable with the range of results found up to now by different groups and 
which lie between 133° and 160° ILplME^gl . 



2. First neighbour distances, pair correlation functions and structure factor 

The second significant effect mentioned at the beginning of this section and which shows 
up as soon as the CP simulation starts, concerns the evolution of the mean distances between 
the Si atoms and their oxygen nearest neighbours. Hence, once the CP dynamics starts, we 
have noticed, by comparison with the BKS results, an elongation of the mean distances be- 
tween the Si atoms and their BO nearest neighbours and a shortening of the mean distances 
between the Si atoms and their NBO nearest neighbours (see Fig. 0). As can be seen in 
Fig. |], these mean distances for the BKS simulated glass have identical values (~ 1.62 A) 
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while there exists an important splitting in the CP case which takes place after a very short 
relaxation (see the inset of the same figure) and remains remarkably stable during the whole 
CP trajectory (the CP Si-BO and Si-NBO mean distances are equal to 1.65 A and 1.58 A 
respectively). 

The shortening of the Si-NBO distance and the lengthening of the Si-BO distance are 
also reflected by the Si-BO and Si-NBO pair correlation functions plotted in Fig. [5[ These 
results are in very good agreement with MD [I2| and RMC studies reported for the 
NS4 glass as well as with the experimental and MD data reported for other crystalline and 
amorphous sodo-silicates n|,3G[]. One has to remark that, in the quoted MD simulations 



fl2| , [14|| the authors have used three body potentials while the potential used in the present 
BKS simulation is a pair potential one. This may be the reason why the quoted MD 
calculations exhibited different values for the Si-BO and Si-NBO mean distances. 

We have also investigated what happens with the mean distance between a NBO atom 
and its Na nearest neighbour. As shown in Fig. || it increases once the CP simulation 
begins. We recall that two Na atoms have been identified in the nearest neighbourhood 
of each NBO atom and the values given in Fig. || are the averages of the two NBO-Na 
distances. The lengthening of this distance which occurs also very shortly after the plugging 
in of the CP description reflects, together with the effects described above, the ability of 
the CP description to change the local environment around the oxygen atoms even at low 
temperature. 

The pair correlation functions (PCF) g a ,p{. r ) (®, f3 = Si, O, Na) and the corresponding 
integrated coordination numbers (CN) are shown in Fig. |7|for both BKS and CP simulations. 
We note that the addition of the sodium oxide does not change the positions of the peaks 
involving the network forming atoms (Si-Si, Si-0 and O-O, Fig. 0(a), (b) and (c)) which 
remain almost the same as those found for pure silica and the two descriptions give overall 
similar results. Hence the first peaks in the Si-0 correlation function exhibit maxima at 
values of r close to 1.62 A, the corresponding experimental value calculated from EXAFS 



spectra |47j being 1.617 A. The Si-0 CN being equal to 4, there is no doubt that the Si04 



tetrahedron remains the basic unit in spite of the presence of the alkali oxide. 

When we compare the PCF involving the modifier cation, we note some differences 
between the results of the BKS and the CP simulations. Firstly, from the analysis of the 
Na-0 PCF (Fig. 0(e)), it follows that the BKS Na-0 first peak is located at « 2.2 A while 
the CP one at ~ 2.3 A. This shift of the average first neighbour Na-0 distance towards 
higher values of r in the CP description is of course consistent with the increase observed 
in the time evolution of the NBO-Na mean distance (see Fig. ||). The CP average first 
neighbour Na-0 distance shows a very good agreement with previous XAFS measurements 
26| , |2"7j| which gave an average value of 2.32 A. 

Concerning the Si-Na PCF (Fig. 0(d)), we can notice that its shape is not similar for the 
classical and ab initio descriptions. After an almost identical first peak located at m 3.4 A, 
the BKS function presents a kind of plateau for r varying from « 4 to 4.8 A. During the CP 
simulation, one can see that this plateau disappears. Since this plateau is more pronounced 
in one of the two BKS samples, it could simply be an artefact of the unrealistic quench 
rates that one has to impose usually in MD simulations. On the other hand, this plateau 
has been observed also in a classical MD simulation using the same BKS potential for one 
NS4 sample containing 648 atoms H3] which shows that it is not due to finite size effects. 
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The first peak position obtained here is in agreement with the RMC |24J] fits and previous 
classical MD results but it is lower than that estimated from sodium K-edge XAFS 

experiments (i.e. 3.8 A) p^j27[ . 



From the same XAFS experiments |26],|27[ , the Na-Na average distance have been refined 
at 3.2 A. In the present work, the first peak of the Na-Na PCF has a maximum at m 3 A for 
both BKS and CP case (see Fig. 0(f)). In the BKS case, the Na-Na PCF shows a splitting 
which seems to be smoothed in the CP case. However, it is difficult to extract really useful 
information concerning the Na-Na correlation given the poor statistics which is due to the 
small number of Na atoms in the system. 

Figure |8| shows the total static structure factor computed from the CP and BKS sim- 
ulations compared with the experimental points that we have extracted from a neutron 
diffraction experiment ||24j| . The differences between the curves corresponding to the two 



simulations are indistinguishable and they agree reasonably with experimental data. Nev- 
ertheless it is interesting to notice that the microscopic discrepancies between the two de- 
scriptions, reported above, are totally smoothed out in the total structure factor. Thus, the 
good agreement shown by the S(q) is not sufficient to draw conclusions on the validity of a 
model and a more detailed analysis of the structure is required. 

B. Electronic properties 

In this subsection, we analyze the electronic properties of the two NS4 ab initio samples 
in terms of the electronic density of states and the Hirshfeld charge variations. These 
properties were computed for the final configurations of the CP runs and compared to the 
ones obtained for the initial configurations (corresponding to the final configurations of the 
classical BKS simulations). We compare the corresponding NS4 results with those obtained 
for the ab initio silica glass sample described in Ref. |29| . 



1. Electronic density of states 

The Kohn-Sham (KS) eigenvalues are computed for both NS4 samples and the resulting 
electronic density of states (DOS) is depicted in Fig. |](a) for the initial configurations of 
the CP runs, and in Fig. 0(b) for the final ones. In Fig. 0(b) , the KS densities of states 



for NS4 are compared to the density of states of the silica glass sample studied in Ref. |f29 
However, since the present NS4 study has been carried out using the LDA-B-LYP functional 
for the exchange and correlation energy term, we have also computed the KS eigenvalues of 
the silica glass using the same approximation (these eigenvalues were originally computed 



within the LDA description in Ref. [29]). We have found only very slight differences between 
the Si0 2 DOS using LDA and LDA+B-LYP. At m 335 K, the band gap is found to be equal 
to 5.01 eV with the LDA+B-LYP description, which is much lower than the experimental 
value for amorphous silica 9 eV). This underestimation of the band gap is a well-known 
failure of DFT in the local density approximation |4^ . 



The band gaps obtained for the two NS4 samples at ~ 400 K are equal to 2.77 eV and 
to 2.86 eV, respectively (Fig. |5|(b)). For comparison, the band gaps obtained for these two 
samples in their initial configurations were equal to 2.16 eV and 2.17 eV, respectively (Fig. 
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|9](a)). Hence, the refinement of the NS4 glass structure by the Car-Parrinello dynamics leads 
to an increase of the band gap. By comparing with the value obtained for the silica glass, we 
note that the gap is significantly reduced by the introduction of Na atoms, as had already 
been deduced from UV absorption experiments |50| , |5l|| and suggested in LCAO calculations 



52]. 

The DOS for the two NS4 samples are very similar and show common features with the 
DOS of Si02 (Fig. |5|(b)). They mainly consist in oxygen 2p non bonding orbitals (lone pair) 
for the highest occupied states, in bonding states between silicon sp 3 hybrids and oxygen 2p 
orbitals (from -10 eV to -4 eV for Si02 and from -11 eV to -5.4 eV for NS4) and in oxygen 
2s orbitals for the lower energy band. The empty band is mainly derived from weak anti- 
bonding conduction states. Note that in our NS4 samples, we take explicitly into account 
semi-core electrons for the Na atoms. These electrons give rise to a low energy band at ~ 
-48 eV which has not been considered in Fig. |^. 

Murray et al. computed the DOS for sodium silicate model glasses with various 
compositions, using a LCAO computational method. Since they used atomic orbitals in 
their model, they could assign different contributions in the bands to states belonging to 
specific atoms. For instance, the NBO atoms gave rise to subpeaks or shoulders to the right 
of the main O 2s and O 2p bands. In the Kohn-Sham description, the orbitals are not 
atomic orbitals and consequently such an analysis of the bands is not possible. However, 
from the comparison of our DOS and the DOS reported in Ref. |52|, we note that they are 
surprisingly similar and therefore, we can attribute the peak just below the Fermi level and 
the peak at ~ -16.5 eV to the NBO atoms. In Ref. []53| |, the reported value for the chemical 
shift between the BO and the NBO ls-photoelectron bands measured by XPS experiments 
was about 2.45 eV. We find values of ~ 2.5 eV for the O 2s bands and of ~ 2.0 eV for the 
O 2p bands. 

By comparison, the bands for the initial configurations - corresponding to the structure 
given by the BKS potential (Fig. |S|(a)) - show only very slight differences with the bands 
computed at the end of the CP runs. The O 2s bands are slightly modified and the peak 
just below the Fermi level attributed to the NBO atoms seems to be shifted towards lower 
values after the CP structure refinement. This is consistent with the Si-NBO bond length 
shortening that occurs after the CP is switched on. 



2. Atomic charges 

It is well-known that the charge of an atom can never been defined uniquely and it is not 
subject to experimental measurement. Nevertheless, it is possible to use approaches such as 
Mulliken population analysis |54] or Hirshfeld determination of the density deformation [55 



in order to discuss bonding or to compare the electron distributions in different systems using 
the same description. Here we have chosen to use the Hirshfeld description to determine 
variations of the atomic charges between the NS4 samples and the silica samples. In this 
scheme, the total electronic charge of a bonded atom is given by: 



Qi = - J 5pi(r)dv 



where Sp(r) is the atomic deformation density, defined as the difference between what Hirsh- 
feld calls the "charge density of the bonded atom" and the "atomic density", i.e. the density 
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difference between the free atom and the bonded atom (see Ref. |55[ for a more complete 
definition). Adding the charge of the nucleus Z^ gives the net atomic charge: 



qi = Qi + Zi. 

This description, based on the integration of the atomic deformation densities, has been used 
to compute atomic charges that appear reasonable in magnitude and show variations 



consistent with the accepted ideas about electronegativity differences between atoms and 
groups. However, the absolute values obtained using this approach tend to be smaller than 
values obtained by other methods, such as Mulliken analysis, and only variations of these 
charges have been examined here. The calculations have been performed at K with fully 
optimized structures in the LDA+B-LYP description and the results are presented in Table 

I 

By comparing the Hirshfeld net charges obtained in the two NS4 samples, we first notice 
that the values are extremely similar, the differences being within the error bars. We have 
also compared these charges with the ones computed on the initial configuration of the two 
CP runs and we have observed no significant differences, given the large error bars. 

In the NS4 glasses, the charges on the NBO atoms are, on average, more negative than 
the ones on the BO atoms, which is a generally admitted result. Moreover we note that 
the differences between the charges on the silicon atoms depend on the character of the 
tetrahedra (Q n ). The positive charge on the silicon atoms decreases as n decreases, i.e. the 
proximity of the NBO and of the Na atoms has a direct effect on the silicon atom and thus on 
the tetrahedra. This effect is known as a "non-localized effect" of the Na atom on the silica 
network |51|j53| , |56[| and has already been shown in ab initio molecular orbitals calculations 



on sodium silicate clusters |5T||56| • But even if the electronic properties of the clusters were 



in agreement with experimental results from UV absorption and XPS measurements [p0| , |53 
some of the geometrical features were not correctly reproduced (in particular NBO-Na bond 
lengths). In our simulation, this "non-localized effect" gives rise to the increase of the Si-BO 
bond length and to the decrease of the Si-NBO one (as presented in subsection |3 A|) , while 
the Na-NBO bond length is in agreement with experimental values. We also observe that 
the increase of the Si-BO bond length is more pronounced in Q3 and Q2 tetrahedra than in 
Q4 tetrahedra, in parallel with the increase of the Si net charges. 

On the other hand, comparison of the average charges found on the BO atoms between 
the NS4 samples and the Si02 sample shows that the BO atoms in Si02 have a more negative 
charge than the BO atoms in NS4. As a direct consequence, the Si atoms in Si02 are less 
charged than the corresponding Si atoms in the Q4 conformation in the NS4 samples. This 
effect is also very probably a consequence of the highly non- localized effect of the Na atoms 
on the silica network. 



4. CONCLUSION 

We have presented a structural and electronic analysis of a sodium tetrasilicate (NS4) 
glass by using a combination of classical and ab initio molecular dynamics simulations. 
The electronic structure has been computed in the framework of the Kohn-Sham density 
functional theory within the generalized gradient approximation using a B-LYP functional. 
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To our knowledge, this is the first study of a binary silicate glass carried out in the frame of 
first-principles molecular dynamics calculations. 

The scheme employed for the preparation of the ab initio NS4 samples economizes a con- 
siderable amount of computation time as was already shown in recent work on the structural 
and electronic properties of a silica glass P5| . As in the silica study, the CP dynamics has 
presented a remarkable stability. In order to validate this approach, we have performed it 
twice and we have obtained the same behavior both times. 

The structural features of the NS4 samples have been studied in terms of pair correlation 
functions, bond angle distributions, Q-species, bridging to non-bridging oxygen ratio and 
structure factor. We note that, once the CP dynamics was plugged in, some important 
structural changes occurred. Hence the shortening of the Si-NBO distance, the lengthening of 
the Si-BO and Na-NBO distances appearing immediately after the start of the CP simulation 
are in perfect agreement with the experimental data. These results validate our preparation 
method and show the ability of the CP description to refine the local environment of the 
atoms even at low temperature. 

The electronic structure of the ab initio samples thus obtained has been analyzed in 
terms of the density of states and the atomic charge variations. The results show that the 
introduction of sodium atoms into the silica network lowers the band gap and that it has a 
highly non-localized effect on the charges of the atoms in the network. 

In conclusion, we have generated - at low computational cost - a fully ab initio sodium 
tetrasilicate glass which shows structural and electronic characteristics in very good agree- 
ment with experimental results. This work is a first step towards accurate studies of detailed 
electronic and vibrational properties of the sodium tetrasilicate glass. 
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FIGURES 
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FIG. 1. Ionic temperature evolution during the last 5 ps of the classical MD simulation (dotted 
line) and the full Car-Parrinello MD simulation (solid line). In the inset, a zoom of the beginning 
of the Car-Parrinello MD simulation is depicted. 
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FIG. 2. Time evolution of the mean O-Si-0 (a), Si-O-Na (b) and Si-O-Si (c) bond angles 
during the last 5 ps of the classical MD simulation (dotted lines) and the full Car-Parrinello MD 
simulation (solid lines). 
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FIG. 3. The O-Si-0 (a), Si-O-Na (b) and Si-O-Si (c) time-averaged angle distributions from 
the Car-Parrinello MD simulation (solid line) and from the classical MD simulation (long dashed 
line), respectively. 
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FIG. 4. Time evolution of the Si-BO and Si-NBO mean distances during the last 5 ps of the 
classical MD simulation and the full Car-Parrinello MD simulation. 
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FIG. 5. First peaks of the Si-BO (solid line) and Si-NBO (dashed line) pair correlation 
functions of the NS4 glass computed from the classical MD (a) and the Car-Parrinello MD (b) 
simulations. 
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FIG. 6. Time evolution of the NBO-Na mean distances during the last 5 ps of the classical 
MD simulation (dashed line) and the full Car-Parrinello MD simulation (solid line). 
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FIG. 8. Comparison between the experimental neutron diffraction structure factor (the cir- 
cles) that we have extracted from reference fl24j] and the computed static structure factors from 
Car-Parrinello MD (solid line) and classical MD (dotted line) simulations. For the calculations of 
the structure factors, the scattering lengths b$i = 4.149 A , bo = 5.803 A and b^a = 3.63 A were 
used. 
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FIG. 9. Upper graph: Kohn-Sham electronic density of states of the initial CP configuration 
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TABLES 





Na 2 0-4Si0 2 






sample 1 


sample 2 


sample from Ref 29| 


BO 


-0.089 ± 0.011 


-0.089 ± 0.009 


-0.109 ± 0.007 


NBO 


-0.245 ± 0.020 


-0.249 ± 0.015 




Na 


0.076 ± 0.038 


0.082 ± 0.034 




Si (all atoms) 


0.240 ± 0.035 


0.240 ± 0.041 


0.218 ± 0.010 


Si (Qi) 


0.262 ± 0.017 


0.270 ± 0.012 


0.218 ± 0.010 


Si (Q 3 ) 


0.222 ± 0.012 


0.212 ± 0.022 




Si (Q 2 ) 


0.125 


0.153 ± 0.005 





TABLE I. Average Hirshfeld atomic net charges computed from the fully optimized structures 



of the two NS4 samples and of the Si0 2 sample of Ref. [29| at K 
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